
!-------------------------------------------------!
!                                                 !
! This program calculates the transmission        !
! function, current-potential curve and projected !
! densities of states for a diatomic molecule     !
! connected to two monoatomic wires               !
!                                                 !
! This module calculates the density matrix by    !
! the combination of a complex integration and    !
! an integration over the real energy axis        !
!                                                 !
!-------------------------------------------------!
! Written by Cristian G. Sanchez                  !
! August 2007                                     !
!-------------------------------------------------!

!-------------------------------------------------!
! The complex path used for the integration of    !
! Gdev to produce the density matrix is a circle  !
! that crosses the real axis at mu_right, the     !
! of the two electron chemical potentials, and    !
! at int_bottom, which must be lower that the     !
! lowest band edge.                               !
!-------------------------------------------------!

module integration

use control
use matrixmod
use gaussquad
use gdevice

implicit none

real(dp), private, allocatable :: quad_points_path(:)  ! evaluation points of quadrature rule
real(dp), private, allocatable :: quad_weights_path(:) ! weights of quadrature points
real(dp), private, allocatable :: quad_points_axis(:)  ! evaluation points for quadrature in real axis
real(dp), private, allocatable :: quad_weights_axis(:) ! weights for axis quadrature

real(dp), private, allocatable :: quad_points_driver(:) ! evaluation points for driver integration
real(dp), private, allocatable :: quad_weights_driver(:) ! weights for driver integration


!real(dp), public ::
integer, public :: algo1

public :: path_integrator
public :: initialize_integration
public :: finalize_integration

 contains

  subroutine initialize_integration()
    real(dp), allocatable :: work(:)
    integer :: info,i

	
		
    !write(*,*)'Initializing path integrator'
    allocate(quad_points_path(path_steps))
    allocate(quad_weights_path(path_steps))
    allocate(work(path_steps))
    call gauss_rule(0,path_steps,quad_points_path,quad_weights_path,work,0.0_dp,0.0_dp,'N',info)
    if (info /= 0) stop 'Gaussian quadrature initialization failed'	
    do i=1,path_steps
      quad_points_path(i) = (quad_points_path(i)+1.0_dp)*pi
    enddo
    deallocate(work)
	
    ! initializing real axis integrator
    allocate(quad_points_axis(real_steps))
    allocate(quad_weights_axis(real_steps))
    allocate(work(real_steps))
    call gauss_rule(0,real_steps,quad_points_axis,quad_weights_axis,work,0.0_dp,0.0_dp,'N',info)
    if (info /= 0) stop 'Gaussian quadrature initialization failed'	
    do i=1,path_steps
      quad_points_axis(i) = (quad_points_axis(i)+1.0_dp)*(mu_left-mu_right)/2.0_dp+mu_right
    enddo
    deallocate(work)
		
	
  end subroutine 

  subroutine finalize_integration()
    deallocate(quad_points_path)
    deallocate(quad_weights_path)
    deallocate(quad_points_axis)
    deallocate(quad_weights_axis)
  end subroutine finalize_integration

  complex(dp) function path_integrator()
    integer      :: i
    real(dp)     :: theta,r,x0,y0,dx,dy,energ
    complex(dp)  :: z, dz
    complex(dp) :: f(N,N),fint(N,N), iimag
    dimension :: path_integrator(N,N)

    iimag=cmplx(0.0_dp,1.0_dp)	
	
    path_integrator = 0.0_dp
	
    r  = (mu_right - int_bottom)/2.0_dp
    x0 = (mu_right + int_bottom)/2.0_dp
    y0 = 0.0_dp
    fint = 0.0_dp

  do i=1,path_steps
		
    theta = quad_points_path(i)
    z  = cmplx(x0+r*cos(theta),y0+r*sin(theta),dp)
    dx = -r * sin(theta)
    dy =  r * cos(theta)
    dz = dx + iimag * dy
    f  = G_dev(z)
    fint = f * dz
		
    path_integrator = path_integrator + fint*quad_weights_path(i)

  enddo

  path_integrator = path_integrator * pi /(2.0_dp*pi*iimag)

  fint = 0.0_dp

  if (abs(mu_right-mu_left)>=1.0D-12) then
    do i=1,real_steps
      energ = quad_points_axis(i)
	fint = fint + Dleft(energ)*quad_weights_axis(i)
   enddo
      fint = fint * (mu_left-mu_right)/2.0_dp
      path_integrator = path_integrator + fint
  endif
		
	
  end function path_integrator

end module integration














